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ABSTRACT 


Solutions  to  the  Resistively  Shunted  Junction  (RSJ)  Equation: 

0 4>  + ( 1 + ycos  <J))  <{>  + sin  <J>  = a + km  sin  wt 

for  small  8 ( = .001)  have  been  approximated  using  the  SDRIVE  integration 
package.  Various  graphic  displays  are  used  to  examine  the  output,  including 
plots  of  <j> , di/dt,  and  sin$  as  functions  of  time;  Poincare  diagrams;  and  plots 
of  the  Lienard  coordinate 

z = 8 <j>  + <{>  + y sin  (j> 

which  has  a close  connection  with  the  "slow  manifold",  as  a function  of  $ , of 
sin  <j> , and  of  time.  Integration  is  performed  by  separation  of  the  second-order 
equation  into  a coupled  pair  of  first-order  equations,  and  numerically 
integrating  with  respect  to  time. 

Several  cases  have  been  examined,  for  y < 1,  representing  quiet  behavior  of 
resistively  shunted  thermometer  oscillator  devices.  The  report  is  an  archive 
record  of  program-test  data. 

A case  of  "jump"(  voltage-spike  ) oscillator  performance  for  y = 1.5  has 
been  simulated  in  considerable  detail,  principally  as  a test  of  the  integrator. 
Parallel,  independent  computer  results  (Sanders  and  van  Veldhuizen,  Amsterdam) 
were  available  for  comparison.  This  case  is  of  considerable  mathematical 
interest,  and  values  of  | y | >1  may  also  occur  in  the  RSQUID  thermometer. 
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INTRODUCTION 


We  are  interested  in  numerical  approximations  to  the  solutions  of  the 
non-linear  equation 


0 <j)  + (1+y  cos<f>)<j>  + sin  <{>  = a + kco  sin  ait  , 

which  is  a model  for  simulating  the  behavior  of  the  circuit  consisting  of  a 
series  connection  of  an  inductor,  a resistor  and  a Josephson  junction  (used  as  a 
noise  thermometer).  Here  the  variable  <\>  is  the  time-dependent 
quantum-mechanical  phase  difference  across  the  Josephson  barrier,  and  the 
parameters  represent  combinations  of  circuit  parameters,  as  well  as  the 
amplitude  and  frequency  of  the  applied  signal.  We  concentrate  on  the  conditions 
of  small  0 as  is  the  case  for  the  noise  thermometer. 

The  value  of  the  parameter  y is  crucial  to  the  behavior  of  solutions.  If 
| Y | <1 , we  have  the  so-called  "quiet"  case.  Taking  |y|>1  gives  "jump" 
behavior.  We  will  examine  this  dichotomy  in  more  detail  later,  but  for  the 

present  consider  the  following:  If  |y|<l,  the  coefficient  of  <j>  is  always 

positive  and,  hence,  acts  as  a damping  force.  However,  for  |y|>l  this  term 

1_ 

changes  sign  (at  cos“l  (“y))  and  switches  from  damping  to  forcing.  This 
causes  rapid  "jumps"  in  the  wave  form.  Physically,  the  system  develops  short, 
high-amplitude  pulses  in  its  voltage  output. 

For  the  quiet  case  |y|<l,0  small  and  a>l  , intuition  about  the  general 
solution  can  be  gained  by  looking  at  the  reduced  first-order  autonomous 
equation 


(1+Y  cos<j> ) ^ = a - sincj) 

<j>(o)  = <j30 

This  is  integrated  directly  for  t as  a function  of  £ , giving  a periodic 
motion  with  period 


2tt 


/ a ^-l 

The  next  step  is  to  generalize  to  the  first  order  driven  equation  for  a)  > 1 
(1+Y  coscj>  )<p  + sin  $ = a + kco  sin  cot 


Using  non-linear  variation  of  parameters,  averaging,  and  letting  y -*■  0, 
Sanders  I11]  shows  that  there  is  a (winding)  solution  with  period 


2tt 

T = 

/ az  - J0Z  (k) 


so  long  as 


ct 

JqOO 


>1 


Here  J0  is  the  zero  order  Bessel  function.  The  Bessel 
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function  behavior  is  borne  out  by  experimental  observation.  Our  numerical 
simulations  show  that  the  above  is  a good  estimate  of  the  period  in  the  second 
order  case,  so  long  as  3 is  small.  The  simulation  also  gives  additional 
detailed  information  about  the  wave  form,  which  is  quite  smooth  for  |y|<|. 

In  case  | y | >1  the  situation  is  very  different.  In  order  to  see  why,  first 
notice  that  defining 

z = 3b  + ( 4>  + Y sincfO 
gives  the  Lienard  co-ordinate  formulation 

z = a - sinb  + km  sin  wt 

3b=  z -(b+Y  sine}))  . 

The  curve  z = <j>  + ysin<j>  in  the  <f>  - z plane  is  called  the  slow  manifold,  [33] 
because  <j>  = 0 on  this  curve.  If  we  take  $ modulo  2tt  , the  curve  has  critical 
points  at  the  two  values  of  cos--*-  ("7)  and  an  inflection  point  at  <p  = tt 


4> 

Suppose  now  that  a>l,  and  that  we  choose  initial  values  of  <j>  and  z on  the  slow 
manifold  (say  at  (0,0)).  The  "average"  value  of  z over  one  cycle  is  positive, 
so  that  z tends  to  increase.  However,  for  points  (b,z)  above  the  slow 

• • 

manifold,  6 is  positive  while  for  points  below  <j>  is  negative.  Hence, 

solution  trajectories  tend  to  be  attracted  to  the  slow  manifold  and  move 
slowly  up  it  until  the  first  critical  point  is  reached.  On  this  part  of  the 
trajectory  d<p/dz  behaves  approximately  like 

(z  ~ (?+Y  sin)  )/S 
" • 
z 3b  + (1+7  cost?)b 

z ~ (b+Y  sin  ?) 

32b  + 3( 1+7  cos)b 

Since  3 is  small,  the  denominator  is  close  to  zero,  and  for  (z,p)  above  the 
slow  manifold  db/dz  becomes  large.  The  solution  "jumps"  across  to  the  other 
branch  of  the  slow  manifold.  (This  description  can  be  made  rigorous).  The 
form  which  the  jump  takes  depends  very  strongly  on  the  initial  conditions.  It 
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is  even  possible  for  the  solution  trajectory  to  reach  the  critical  point  and 
loop  back  to  the  first  branch  before  jumping  across  to  the  second  branch. 


The  time-dependent  behavior  of  the  wave-form  in  this  situation  can  be  seen  in 
Figs.  23-37  below. 

Using  an  asymptotic  expansion,  Sanders  t11]  obtains  expressions  for  the  time 
at  which  the  jump  occurs.  If  k=0  (the  autonomous  second  order  equation),  and 
y is  close  to  1 , the  period  of  the  solution  is  approximated  by  T = To  + AT 


where 


2u 


To  = 


/ a^-1 


and 


AT  = 2 
2 


(1-y)' 


In  two  driven  cases  which  we  have  examined,  with  larger  values  of  y-1  (y=  0.5 
and  y = 1.5),  this  formula  yields  too  large  a value  for  the  observed  AT;  the 
time-averaged  limit  cycle  is  better  approximated  by  the  Bessel-function 
expression  mentioned  above. 
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HISTORICAL  SYNOPSIS 


^ 1 
The  second-order  RSJ  equation,  in  autonomous  form,  was  proposed  by  McCumber  [x] 

and  Stewart  [2]  on  the  basis  of  the  Josephson  relations: 

V (t)  = 34>  /3 1-  h/2e  Isc(t)  * Ic  sin  <t> 

combined  with  Kirchoff's  laws  for  the  circuit  in  which  the  tunneling  element  is 
connected.  (Here  V and  Isc  are,  respectively,  the  voltage  across  and  the 
supercurrent  through  the  tunneling  barrier;  <j>  is  the  superconducting-phase 
difference;  and  Ic,  the  critical  current).  Their  objective  was  to  explain  the 
(time-averaged)  dc  current-voltage  characteristic  curve  which  they  observed  for 
the  circuit. 

Kamper  [3],  and  Kamper  and  Zimmerman  [4]  have  suggested  that  a point  contact 
(the  Josephson  junction)  facing  an  inductive  loop  in  series  with  a very  low 
resistor  might  serve  as  a unique  form  of  noise  thermometer.  Soulen  and  Giffard 
[5 ] applied  the  RSJ  model  to  describe  the  behavior  of  the  circuit  in  which  the 
Bessel  function  dependence  was  derived.  A recent  paper  by  Gallop  [6  ] surveys 
the  importance  of  Josephson  SQUID  noise  thermometry,  as  an  absolute  temperature 
calibration  standard. 

Peterson,  Soulen,  and  Van  Vechten  [7],[8],[9]  examined  a first-order-accuracy 
solution  to  the  time-dependent  driven  equation.  (A  recent  manuscript  by  Sanders 
t11]  clarifies  the  same  picture  for  second  order.)  This  provides  a partial 
description  of  the  frequency  perturbations  in  the  relaxation  oscillation 
observed  when  the  drive  amplitude  is  varied.  In  addition  to  the  relaxation 
frequency  dependence  on  drive  amplitude,  through  a Bessel  function  expression 
(which  is  expected  from  the  model,  and  observed  experimentally),  a weak  upward 
or  downward  slope  is  found  experimentally.  This  may  be  simulated  partially  by 
expanding  the  Bessel-function  solution  to  the  first-degree  equation,  and 
partially  by  the  inclusion  of  interactions  with  the  (external)  tank  circuit. 
Puzzling  problems  remain  such  as  the  observed  occurrence  of  "doubled”  peaks;  R. 
Soulen,  Private  Comm.,  [ 1 1 A]  and  it  is  clear  that  neglect  of  the  second 
derivative  leaves  the  treatment  incomplete. 

Schlup  (1979)  [10]  has  offered  a relatively  complete  model  for  behavior  of  the 
small-beta  autonomous  case.  The  treatment  by  Sanders  [ 1 1 ] links  the  autonomous 
relaxation  oscillation  to  a Hopf  bifurcation  from  the  solutions  on  the 
zero-voltage  step  (that  is,  when  the  bias-current  term  alpha  has  values  below 
one. ) 

Driven  Case  for  beta  large 

The  case  of  beta  large  (8»  1 ) has  recently  been  rather  broadly  treated. 

Schlup  (1974)  [12]  provided  a numerical  solution  for  the  time  dependence,  in 
this  the  near-Hamiltonian-oscillator  case.  His  paper  was  the  first  to  include 
specific  effects  from  the  y term.  Belykh,  Soerensen,  and  Petersen  I and  II 
[i3,i4]  have  provided  a survey  of  possible  solutions  for  the  autonomous  and 
driven  cases,  respectively. 


*The  abbreviation  RSJ  (Resistively  Shunted  Junction)  equation  is  in  common  use 
for  this  equation,  which  may  variously  appear  without  the  3 4)  , y cos  <J>  , or 

k a)  sin  ojt  terms.  The  term  R-SQUID  noise  thermometer  refers  to  the  circuit 
described  above. 
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Connections  of  the  case  y = 0 and  0 > 1 with  bifurcations,  chaos,  and 
oscillator  instability  have  been  examined  by  Packard,  Huberman,  et  al. 
[15,16,1?];  in  detail  by  R.  Kautz,  [18,19,20],  and  by  Pederson  and  Davidson 
[ 2 1 ] , who  looked  at  amplifier  stability  over  a broad  range  of  frequency  and 
amplitude  parameters.  This  case  has  important  applications  to  phase-locked 
f requency-to-dc-voltage  converters,  and  to  microwave  parametric  amplifiers. 

Driven  case  for  beta  small,  and  gamma  non-zero. 

The  detailed  parametric  treatment  of  the  RSQUID*  noise  thermometer  has  been 
begun  by  Soulen,  Park,  Seppa,  and  Van  Vechten  [22,23,24]  . Careful  estimates 
of  realistic  parameter  ranges  have  been  made  by  Van  Vechten  [2d]. 

Extremely  fast  transitions  between  parts  of  the  slow  manifold  appear  at  once 
if  gamma  is  assumed  larger  than  one.  These  render  the  solution  of  the  first 
degree  equation  at  best  an  unsatisfactory  approximation  to  the  full  wave  form. 
Cushman  [26 ] has  completed  a global  treatment  of  existence  of  periodic 
solutions  for  the  second-order  autonomous  equation,  based  upon  a general  model 
of  flows  in  phase  space.  Important  results  from  this  case  include  indications 
of  the  strong  global  stability  of  solutions  in  the  autonomous  case,  for  all 
values  of  gamma;  and  rapid  convergence  of  these  to  the  periodic  solution. 

Sanders  [ 1 1 ] has  offered  an  approach  to  the  ’'jump"  ( y > 1 ) cases,  which  is 
inspired  by  the  work  of  Levi  [ 2 8 ] , and  which  proceeds  from  recent  successful 
treatments  of  the  van  der  Pol  equation  by  Grasman  and  collaborators  [29,30]. 
These  techniques,  in  turn,  are  based  upon  Smale's  work  in  differential 
topology. 

On  the  basis  of  this  wide  spectrum  of  interest  in  the  mathematical  features  of 
the  driven  (forced)  equation,  and  its  physical  consequences  as  a simulation, 
we  have  been  led  to  an  approach  which  stresses  detailed,  accurate  numerical 
integration  combined  with  varied  graphics  displays. 


RSQUID  = Resistor  - Superconductive  Quantum  Interference  Device 
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THE  ALGORITHM 


The  numerical  procedure  for  integration  has  been  to  separate  the  second-order 
equation  into  two  coupled  first-order  equations: 


Yl  = Y2  / 8 


(1)  Y2 


- (1  + Y cos  Yl)*Y2  /8  + a - sin  Yl  + kw  sin  wt 


so  that 


Y1  = (j> 


Yl  = 4> 


Y2  = 84>  Y2  = 6<j> 

This  is  not  a unique  separation.  Another  choice  for  Lienard  coordinates  comes 
from  noticing  that 


8c()  + ( 1 + ycos  $)<j>  , is  an  exact  differential. 


In  that  case 
let 


z = 8$  + ( <j>  + -ysin  $ ) 


, so  that 


(2) 


and 


z = a - sin  cf>  kxu  sin  wt 
8£=  z-($+y  sin  $ ) 


Although  the  first  choice  of  Lienard  coordinates  appears  to  involve  more 
multiplications,  this  choice  is  numerically  more  stable  than  the  second 
because  8 is  small  and  z stays  close  to  (<j>+Ysin  <p).  However,  the  second  pair 
has  some  independent  mathematical  interest,  as  has  been  mentioned. 


Integration  is  performed  by  using  either  SDRIV3  in  single  precision  or  DDRIV3 
in  double  precision  These  programs  utilize  implicit  raultistep 

integration  formulas,  either  of  Gear  or  Adams  Moulton,  to  track  the  solution 
components  of  (1)  or  (2)  to  within  certain  specified  local  error  criteria. 

The  nonlinear  algebraic  equations  which  result  can  be  solved  by  any  of  several 
algorithms,  including  Newton's  method,  in  which  case  the  Jacobian  of  the 
system  can  be  given  explicitly  or  computed  numerically.  Various  output 
options  are  also  available. 


STRUCTURE  OF  THE  PROGRAM 

The  main  program,  CU SHEQ/ DOUBLE , simulates  the  operation  of  the  physical 
system  by  querying  for  parameters  and  initial  value  settings,  through  the 
subroutine  OQUERY / DOUBLE ; driving  the  DDRIVE  integrator  with  its  attached  FI 
(Non-autonomous  terms)  and  JACOBN  (Jacobian)  subroutines;  returning  values  at 
time  intervals  of  S • 2tt  where  S is  an  input  real  number;  and  storing 

w 

the  values  of  t,  Yl  , Y2  in  arrays. 
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These  values  are  written  to  a file  by  the  subroutine  OPRINT /DOUBLE  and  may 
later  be  graphed  in  a variety  of  ways  through  the  program  SPREAD / DOUBLE , 
which  offers  nine  options  of  two-dimensional  plotting.  SPREAD  collects  the 

values  of  <j>  and  cj>  , and  evaluates  sin  * and  z (to  single  precision)  in  one 
loop. 

This  functional  structure  of  the  subroutines  is  illustrated  in  Table  I. 

INPUT  INSTRUCTIONS 

The  operator  is  queried  for  input  values  of 
A = a 


B = 0 
C = Y 

D = 5 = 1 (a  parameter  coefficient  of  the  sin  <j>  term,  which  is  kept  in  reserve 
for  later  approximations) 

E = kw 

2e 

W = a)  , the  circular  frequency  of  the  drive  term,  scaled  to  = — RIC  * 

* 

S , a real  number,  which  selects  the  time  intervals  of  At  = S*2tt  at  which 

u> 

numerical  results  will  be  returned  from  the  integrator;  and  for  initial 
choices  of 

Y1  = <J> 

Y2  = J 
T = t 


TOL  = the  numerical  tolerance  for  the  completed  integration.  Finally,  the 
number  N (=  number  of  points  to  be  returned)  is  requested. 


* See  Appendix  B for  the  circuit 


SUBROUTINE  STRUCTURE  OF  THE  PROGRAM  CUSHEQ 


OPTIONS  (9) 


GRAPHICS  OUTPUT  CHOICES 


The  graphics  output  subroutines  can  be  combined  with  appropriate  choices  of  S to 
produce  a very  large  variety  of  displays.  These  will  only  be  briefly  described; 
they  are  best  illustrated  through  examples  of  use. 

Option  (1):  The  Poincare  map  routine,  of  <Jj  vs  mod  2ir  , is  conventionally  used 

in  non-linear  analysis  to  distinguish  fixed  points,  stable  and  non-stable 
orbits,  etc.  In  the  present  case,  selection  of  S permits  the  time  interval  to  be 
taken  as  the  period  of  the  drive  oscillation  (S  = 1),  or  any  multiple  or 
fraction  thereof. 

Option  (2):  As  mentioned  on  p.  9 , z is  constructed  by  SPREAD  from  the  stored 
values  of  and  Y2.  The  plot  of  <j>  vs  z is  useful  in  distinguishing  transition 

phenomena  for  the  high-speed  ''jump''  intervals  which  occur  for  the  choices  of 
gamma  larger  than  1 , where  the  term  y cos  <j>  can  have  negative  values  greater  than 

1 . 

Option  (3):  $ vs  t is  the  actual  time-dependent  solution  of  the  equation, 

simulating  the  barrier-phase-difference  for  the  Josephson  effect.  <J>  is  plotted 
mod  2tt  , since  otherwise  the  secular  term  in  its  solution  would  cause  it  to  "ramp" 
off  the  page. 

4>  vs  t simulates  the  instantaneous  time-dependence  of  the  barrier 
voltage. 

Option  (4):  The  return  map  of  <p  at  intervals  corresponding  to  the  drive 

period  is  a method  of  exploring  both  the  subharmonic  structure  of  the  solutions, 
and  the  jumps  between  branches  of  the  slow  manifold  which  the  system  may  exhibit. 

Option  (5):  sin  <jj  vs  t simulates  the  instantaneous  flow  of  supercurrent  through 
the  barrier  of  the  Josephson  junction.  There  are  a number  of  interesting 
insights  and  visual  simplifications  introduced  by  the  cylinder  projection  of 
viewing  which  the  sine  function  provides.  The  physical  quantity  <j>  is  not  an 
observable;  it  represents  the  argument  of  a trignometric  function  which  is  an 
observable,  and  $ can,  therefore,  be  thought  of  as  confined  to  motion  on  a 
cylinder. 

Option  (6):  sin  ^ vs  z and 

Option  (7):  z vs  sin  $ illustrate  the  same  effects  as  Option  (2),  but  in 

cylindrical  projection,  relating  z to  the  (observable)  supercurrent. 

• 

Option  (8):  z vs  t collects  the  values  of  <jj  and  <j>  as  returned  by  the  main 
integration,  and  calculates  z , displaying  it  as  a function  of  time.  The 

smoothness  of  this  function,  when  its  individual  terms  are  undergoing  rapid 
derivative  changes  with  repect  to  time,  is  a measure  of  the  performance  of  the 
program  as  a whole. 

Option  (9):  ip/t  vs  t provides  a mechanism  for  averaging  the  advance  rate  of 

the  quantum-mechanical  phase  difference,  a measure  of  the  relaxation,  or  limit 
frequency 

<q>  = lim  _1  ,p(t)dt 

y-HX  x 0 

It  is  calculated  by  computing  Yl/t  values  from  the  output  array,  while 
suppressing  the  mod  2tt  reduction  used  in  plotting  $ , in  option  (3)  above. 
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SCALING  OF  THE  EQUATION  TO  THE  PHYSICAL  PROBLEM 


Time  Scaling 


Drive  and  relaxation  frequencies  are  scaled  to  the  real  physical  parameters  of 
the  thermometer  oscillator  circuit  through  the  "equivalent  plasma  frequency" 


2e 

2tt  ujr  = RIr 


in  circular  measure  , or 


v 


c 


in  direct  frequency  measure  (Hz) 


2e 

The  quantity  is  known  extremely  accurately,  from  voltage  stabilization 

h 

experiments,  to  have  a value  of  483.6  MHz/microvolt  *. 

5 6 

and  for  R = 10“  ohms  (the  shunt  resistor-see  Appendix  B)  with  Ic  = 10”  amps 


vc  = (483.6  x 10  12  ) x 10  ”5  x 10  “6  Hz 

vc  = 4.836  x 103  Hz  or  4.836  kHz  , specifically  when  the 
product  RIC  is  equal  to  10-11  volts. 


Selected  Drive  Frequencies 

A value  for  co  of  200,  scaled  to  the  circuit  parameters  above, 
corresponds  to 


200  x 4.836  x 

10' J 

= 967 

kHz 

Or,  a value  for  u of  800  corresponds 

to 

3.869 

Mhz 

And,  a value  for  oj  of  45.68  is  equivalent  to 

45.68  x 4.836  x 103  = 220.9  kHz 

These  are  the  values  of  drive  frequency  used  in  the  simulations  reported 
below. 


* From  high-precision  experiments  against  the  U.S.  standard  volt  cells,  this 
quantity  is  known  at  better  than  one  part  per  million  accuracy.  The  rounded 
version  to  four  decimal  figures  is  used  here,  for  computational  simplicity. 
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Time-dependent  quantities 


Time-dependent  quantities  in  the  RSJ  equation  are  expressed  in  voltage  units 
for  simulation  purposes: 

thus  a = 1 for  R Ic  = 10“  ^ volts-dc,  in  the  examples  used,  and  therefore 
kw  = 1 corresponds  to  an  rf  input  signal  at  10“^  volts-rf  amplitude  (not  rms). 

Output  <J>  = Y2/3  = 1 has  the  scale  (10“®  volts)  and  sin  <}>  = 1 corresponds  to 
I s c = Ic  = 10"6  amps. 

The  parameter  3 

7 

Soulen  [ ] estimates  a value  of  the  second-derivative  coefficient  3 at  .001,  as 
reasonable  for  point  contact  RSQUID  devices.  We  note  that  a capacitative  value 
of  3 as  small  as  .01  has  been  suggested  for  experiments  on  tunnel  junctions,  by- 
Fulton  [34].  This  does  not  provide  strong  guidance  for  the  case  of  point 
contacts,  however,  because  of  the  different  physical  barrier-layer  present.  The 
value  .001  is  consistent  with  the  curves  of  Fig  2 and  Fig.  5 in  McCumber  ]. 
We  select  beta  positive  to  assure  that  the  second  derivative  term  will  have  a 
limiting  effect  during  the  " jump"  intervals  for  cases  of  y greater  than  one. 

Beta  vs  omega  plane 

If  we  cons^er  the  plot  of  beta  versus  reduced  frequency  u>/wc  which  is  found 
in  Kautz  [ ] as  Fig.  2 of  that  article,  the  value  3 = .001  in  conjunction 

with  the  line  where 

wL  = 1/  wC  , its  resulting  value 

3 “ 1/2  = 31.62 

would  lie  along  the  lower  edge  of  region  VI  extended  to  the  upper  left, 
in  Kautz* s log-log  beta-omega 

plane.  Since  the  values  chosen  for  uj/  wc  , namely  200,  800,  and  45.68,  are 
above  this,  our  simulations  lie  up  in  the  region  between  that  lower-edge  value 
and  the  line  where  R = 1/uC 

3-1  = 1000 

This  is  equivalent  to  saying  that  our  results  have  been  limited  to  a frequency 
region  above  the  junction/circuit  "equivalent"  plasma  frequency.  Such 
inf ormation  is  helpful  when  comparing  behavior  between  different  families  of 
Josephson  oscillators. 
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Choice  of  y 


y from  theory  is  proportional  to  Ic  Le  , where  Ic  is  of  the  order  of  10- 
amps  and  Le,  the  inductance  of  the  coupling  loop  which  the  Josephson  junction 
faces,  is  estimated  to  be  about  2xl0-^  henries. 

If  we  take  the  expression  derived  by  Petersonty]  for  gamma  in  the  RSQUID 
circuit : 

2tt  2e 

Y = — - * 2e  x I^IC  = 2tt  x x Lelc 

= 2-rr  x 483.6  x l012x  2xl0”10  x 1CT6  = 0.6077 

for  the  above  values.  Here  we  see  that  Le>  0.  33  nanohenry  will  suffice  to 
render  y > 1.  y » as  shown  by  Soulen,  Peterson,  and  Van  Vechten  [7 ] can  be 
adjusted  experimentally  to  fall  either  below  or  above  1 . It  is  a parameter 
which  can  be  determined  in  a coarse  manner  within  the  experiment.  Of 
particular  interest  are  the  cases  in  the  transition  region  as  y approaches 
1 from  below.  In  some  cases  there  are  anomalies  in  the  observed  relaxation 
frequency  as  a function  of  drive  amplitude.  Previous  work  has  not  shown 
whether  all  of  these  result  from  an  intrinsic  wave  form  of  the  relaxation 
oscillator,  or  from  its  coupling  to  the  drive  and  detection  circuits. 

A.  Tests  of  Quiet  Cases  with  Smooth  Wave  Forms 

Simulation  of  the  quiet  thermometer  oscillator  has  proceeded  by  examining 
several  cases  for  y < 1 • These  results  can  be  compared  with  existing 
asymptotic  models  (Soulen  [5,7];  Sanders  t11])  and  directly  with  experiment. 
The  graphics  routines  have  been  used  to  simulate  the  waveforms  in  considerable 
detail. 

The  curves  which  follow  in  Figs.  1-10  correspond  to  the  case  y = 0.8  , oj  = 
200,  a = 0.5  , which  corresponds  to  a drive  frequency  of  967  kHz  and  a 
relaxation  frequency  of  about  25  kHz.  Figs.  11-15  are  run  with  the  drive 
frequency  four  times  as  large;  i.e.  3.869  MHz. 


Figs.  1-5  — "Test  30" 

Parameters : 

A = 5,  B = .001,  C = 0.8 ,D  = 1.0,  E = 382,  S = .025,  T = 0,  TOL  = 10-8  , 

W = 200,  3000  points 

This  is  a fairly  low  drive  amplitude,  simulating  a waveform  which  is 
primarily  the  relaxation  oscillation,  modulating  a "carrier"  wave  at  the 
drive  frequency. 


Fig.  1 shows  phi,  the  solution,  as  a function  of  elapsed  time  t.  The 
fine-grain  vertical  oscillations  are  at  the  drive,  or  "carrier"  frequency;  the 
broad  stripes  conform  to  the  relaxation  or  limit  cycle  of  the  Josephson 
oscillator . 
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The  points  on  Figs.  1,  2,  3 have  been  returned  by  the  program  at  the  rate  of 
one  per  1/40  of  the  drive  cycle.  The  numerical  values  have  been  generated  at 
double  precision  (parts  in  108).  Phi  is  plotted  modulo  2tt  , to  prevent  the 
plot  from  ramping  diagonally  off  the  paper. 

Fig.  2 shows  the  time  derivative  of  phi,  simulating  the  instantaneous 
Josephson  voltage  across  the  tunneling  barrier.  Note  that  this  appears  to  be 
primarily  amplitude-modulated  by  the  relaxation  oscillation,  but  that  the 
positive  peaks  in  the  envelope  are  slightly  offset  in  time  from  the  negative 
ones.  The  mean  value  of  phidot  is  expected  to  approximate  A,  or  a ; the 
time  offset  is  a consequence  of  the  value  of  C,  or  y 

Fig.  3 shows  the  sine  of  phi,  simulating  the  instantaneous  super current 
through  the  barrier.  This  is  equivalent  to  taking  the  function  of  Fig.  1 and 
wrapping  it  around  a cylinder,  whose  axis  proceeds  along  the  time  direction. 
The  resulting  "wrapped  ribbon"  is  viewed  from  the  side;  note  that  the  zero  at 
the  bottom  of  Fig.  1 has  now  been  displaced  to  the  center  of  Fig.  3.  Note 
that  the  ribbon  returns  down  the  "far"side  of  the  cylinder  with  steeper  slope 
than  it  rises  on  the  "near"  side.  A time-smoothed  version  of  the  instantaneous 
supercurrent  would  appear  to  be  primarily  amplitude-modulated  at  the 
relaxation  frequency;  Fig.  3 suggests  that  this  would  be  a somewhat  asymetric 
sawtooth  as  a function  of  time,  with  a slow  rise  and  a more  rapid  return  to 
the  negative  limit  of  -1. 

Fig. 4 shows  z versus  phi  for  these  points.  4b  and  4c  illustrate  z vs.  sin  4)  , 
to  suggest  the  cylindrical  projection  which  controls  the  actual  values  of  the 
tunneling  supercurrent.  Fig.  5 shows  the  composite  Poincare  map. 

Each  point  in  these  diagrams  can  be  tracked  to  corresponding  instantaneous 
values  in  Figs.  1,2,3.  These  methods  of  display  are  useful  in  examining 
periodicity  and  stability  properties. 

Figs.  6-10  — "Test  38" 

Parameters : 

Identical  with  "Test  30",  except  E=1402  , about  3.6  times  as  large.  The 
waveforms  become  more  frequency-modulated  in  character. 

Fig.  6 shows  phi  as  a function  of  time.  Note  that  the  fine-grain  vertical 
oscillations  have  been  replaced  by  stripes  moving  along  the  ribbon,  and  that 
some  of  these  stripes  appear  to  cross  over  each  other.  Note  that  the  open 
bands  which  appeared  before  have  now  been  eliminated;  the  ribbon  runs  together 
or  overlaps  on  itself. 

Fig.  7,  corresponding  to  Fig.  2,  now  shows  quite  a different  form  of  voltage 
signal.  The  maximum  amplitude  is  nearly  uniform,  as  in  a frequency-modulated 
waveform.  A sixth  harmonic  of  the  relaxation  oscillation  shows  up.  The 
offset  of  positive-  from  negative-peaks  seen  in  Fig.  2 is  now  reflected  in  the 
interior  "cluster"  detail. 
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Fig.  8.  The  ribbon  winding  of  Fig.  3 has  been  widened  to  the  point  that  it 
overlaps  itself.  At  the  same  time,  the  "carrier"-f requency  detail  is  now 
somewhat  suppressed,  as  we  should  expect  if  the  sidebanding  of  frequency 
modulation  is  superimposed. 

Figs.  9 and  10  correspond  to  Figs.  4 and  5.  Some  changes  of  form  are 
observed. 

Test  of  Higher  Drive  Frequency  and  "Demodulation” 

Figs.  11-15  — "Test  61" 

Parameters : 

A = 5,  B = .001,  C = 0.8,  D = 1.0,  E = 1528  , S = 1.0,  T = 0,  TOL  = 10-4  , 
W=  800,  3000  points 


This  is  roughly  similar  to  the  conditions  of  Figs.  1-5,  except  that  the 
drive  frequency  has  been  increased  by  a factor  of  four,  and  the  amplitude  of 
the  drive  signal  increased  to  produce  roughly  the  same  level  of  amplitude  of 
the  "carrier".  In  addition,  only  one  point  is  returned  per  cycle  of  the 
drive,  instead  of  forty;  thus,  the  drive  signal  disappears  from  view  in  the 
output.  The  horizontal  time  scale  has  been  compressed  by  a factor  of  ten  in 
plotting. 


Fig.  11  shows  phi  as  a function  of  elapsed  time.  The  profile  of  the 
relaxation  oscillation  is  the  dominant  detail.  In  electrical  engineering 
terms,  the  "audio”  signal  has  been  mathematically  "demodulated"  from  the  "rf 
carrier"  and  its  sidebands,  and  is  displayed  here  by  itself. 


Fig.  12  shows  instantaneous  voltage.  Here  we  are  seeing  a sort  of 
stroboscopic  view  of  one  of  the  traces  which  can  be  seen  extending  through 
Fig.  2.  It  gives  only  a limited,  partial  view  of  the  full  waveform  which 
would  be  experienced  if  we  were  seeing  all  phases  of  the  "carrier"  cycles. 


Fig.  13  , a correspondingly  "stroboscopic"  view  of  the  supercurrent  could  be 
correlated  with  one  of  the  lengthwise  threads  along  the  ribbons  in  Fig.  4 or 
Fig.  8.  Note  that  this  is  a continuous  smooth  path  — but  in  Figs.  4 and  8 
there  may  be  crossing  of  these  paths.  Whether  ribbon  twisting  can  take  place 
under  the  proper  conditions  on  C and  E (gamma  and  k omega)  is  an  interesting 
problem  relating  to  the  frequency  stability  of  the  oscillator. 

Figs.  14  and  15  are  smooth  curves,  much  as  expected  for  a stable  system. 
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Tests  at  a Low  Value  of  (a  - 1) 


We  may  select  a = 1.2  in  order  to  emphasize  the  perturbations  on  the 
relaxation  frequency  which  are  caused  by  varying  the  amplitude  of  the  drive 
term.  At  the  same  time,  the  selection  of  oj  = 200  simulates  a drive  frequency 
of  967  kHz,  while  the  chosen  value  of  alpha  corresponds  to  a relaxation 
frequency  of  ~ 5.8  kHz.  We  now  have  quite  a large  band  separation  between 

the  drive,  or  carrier  frequency,  and  the  relaxation  or  limit  cycle  frequency 
of  the  Josephson  oscillator.  Figs.  16-22  below  are  a simulation  of  this  case. 


r 5 8 9 l 

The  theory  by  Peterson,  Soulen,  and  Giffard  *■  > » J , which  is  based  on 
consideration  of  the  first-derivative  term  only,  gives  a value  for  this 
correction: 

J0  2 (k)  1/2 

u)re2_  = a ( 1 - ) x a)c 

a2 

A series  of  curves  has  been  run  for  these  values  of  alpha  and  omega,  over  a 
range  of  k from  1.91  to  8.65,  which  is  designed  to  cover  the  first  three  zeros 
of  the  Bessel  function.  The  relaxation  frequencies,  which  may  be  read  off  to 
better  than  1%  from  the  plots,  do  follow  the  amplitude  dependence  predicted 
from  this  theory  of  the  first-order  equation.  This  corresponds  to  the  period 
of  the  "winding"  solution  of  Sanders,  mentioned  in  the  Introduction. 

Figs.  16-18  — "Test  4" 

Parameters : 

A = 1.2  , B = .001,  C = 0.5,  D = 1.0  , E - 730,  S = 0.20  , T = 0 , TOL  = 10  -4 
W = 200;  3000  points. 

There  are  five  points  returned  per  drive  cycle. 

Here  we  expect  fully  quiet,  smooth  oscillator  behavior.  Because  of  the  low 
value  of  alpha,  different  paths  occur  for  the  off-set  "strobe"  positions 
around  the  drive  cycle.  This  is  particularly  clear  in  the  "wrapped-cylinder" 
paths  of  Fig.  18. 

We  note  parenthetically  that  the  projection  of  Fig.  18  is  taking  the  sine  of  a 
phi-function  which  already  comprises  several  low  harmonics  of  the 
"fundamental"  relaxation  frequency.  "Sine  of  sines"  implies  a Bessel-f unction 
set  of  harmonic  coefficients.  In  this  case,  the  argument  of  the  Bessel 
functions  is  the  drive  amplitude.  Therefore  in  Fig.  18  it  is  not  surprising 
that  we  find  patterns  suggestive  of  very  simple  Fourier-component  mixes,  in 
various  relative  time  phases. 

The  patterns  seen  in  Figs.  17  and  18  bear  close  qualitative  resemblances  to 
the  experimentally  "demodulated"  audio  spectrum  of  an  actual  thermometer,  when 
the  first  six  or  eight  harmonics  of  the  relaxation  frequency  are  allowed  to 
pass  the  audio  filter. 
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Figs.  19-21  — "Test  7” 

Parameters  are  the  same  as  in  Figs.  16-18  above,  except  that  the  drive 
amplitude  is  now  increased  to  1730,  i.e.  2.37  times  as  large. 

At  this  particular  amplitude  the  fundamental  Fourier  coefficient  J^(k)  must  be 
at  a high  value.  Also,  the  paths  overlap  from  one  relaxation  cycle  onto  the 
next,  which  is  only  slightly  the  case  in  test  4 above. 

Fig.  22  illustrates  the  perturbation  of  the  Josephson  relaxation  frequency  by 
the  amplitude  of  the  drive  signal,  which  is  one  of  the  key  problems  in 
physical  variation,  for  which  an  answer  is  sought. 

The  crosses  are  estimated  numerical  values  of  the  relaxation  frequency  sealed 
to  the  critical  frequency  wc  , obtained  from  the  2tt  - crossings  of  phi  vs  t , 
in  the  series  of  which  Figs.  16  and  19  are  examples.  Scatter  of  these  points 
is  about  6 parts  in  10^;  it  comes  partly  from  numerical  roundoff  of  the 
single-precision  integration,  and  partly  from  interpolation  errors  in 
estimating  the  2tt  crossing  points. 

The  circles  are  points  taken  by  interpolating  between  curves  of  Peterson  [8, 
esp.  Fig.  2),  which  is  actually  calculated  for  a = 1.25  , rather  than  u = 1.2, 
the  value  for  which  Figs.  16,  19  ...  have  been  run.  Note  that  the  curves  in 
Peterson's  article  are  plotted  on  a squared  scale;  here  we  show  frequencies 
directly.  The  curve  described  by  the  circles  has  been  normalized  to  the  first 
two  maxima  in  the  crosses  — this  scale  is  shown  at  the  right. 

We  conclude  from  this  comparison  that  the  positions  vs  k of  the  first  two 
maxima  and  minima,  and  the  magnitude  of  the  swings  between,  show  agreement 
within  the  available  limits  of  error.  In  broader  terms,  the  second-order 
simulation  appears  to  be  agreeing  in  detail  with  the  predictions  of  the  first- 
order  model,  which  in  turn  are  known  to  model  a portion  of  the  observed 
experimental  behavior. 

3 . Tests  of  "Jump"  Case(y  ^ 1.5)  with  Discontinuous  Waveform  Derivatives 

The  case  of  gamma  larger  than  one  can  occur  in  R-SQUID  noise  thermometers.  This 
comes  about  by  increasing  the  pressure  on  the  Josephson  point  contact,  so  that 
the  critical  current  Ic  increases,  for  the  same  value  of  exterior  Le  into 
which  the  junction  faces;  gamma  is  proportional  to  the  product  of  Ic  and  Le. 

The  first-derivative  term  in  the  equation  causes  rapid  transitions  when  (1  + y 
cos  $)  changes  sign.  These  rapid  transitions,  on  a greatly  compressed  time 
scale,  give  a very  rigorous  test  of  the  integrator.  While  cos  «■>  is  negative, 

• 

b rises  or  falls  very  quickly  to  an  extreme  value,  then  turns  around  and 
returns  quickly  until  the  solution  is  "caught"  again  as  the  y cos  -j>  term  reaches 
-1  , and  the  system  reverts  to  a "slow  manifold"  solution. 

The  problem  of  "jumps"  of  this  type  has  already  been  studied  for  the  case  of  the 
Van  der  Pol  oscillator  [23,3a]  "Drops"  and  "slices"  occur,  as  presented  by 
Littlewood  [3i,32].  Our  approach  has  been  primarily  through  graphics 
techniques,  using  the  integrator  program  to  simulate  the  physical  evolution  of 
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the  system  as  nearly  as  possible.  We  have  simulated  a description  through  the 
"slow  manifold"  approach,  which  was  suggested  by  Cushman  [33].  We  have  had 
available  to  us  a set  of  independent  calculations  by  Sanders  and  Van  Veldhuizen, 
[27]  using  a different  Shampine-Gordon  integrator,  which  confirm  critical 
details  of  the  program  performance,  especially  for  these  "jump"  details. 

The  parameter  values  are  chosen  as  follows: 

0 = .001  , y = 1.50  , k = 1.0  , a = 11.0  oj  = 45.6800  , ku  = 45.6800 

S has  been  chosen  in  a variety  of  ways,  to  display  different  details. 

This  simulates  a Josephson  thepiometer  with  shunt  resistor  value  R = 10“5  ohms, 
with  critical  current  Ic  = 10”  amps,  driven  at  220.9  KHz  and  at  a drive 
amplitude  k = 1.0,  which  places  its  operation  in  the  "very  weakly  driven"  part 
of  the  Bessel-function  plot  for  relaxation  frequency  (below  the  first  maximum). 
We  note  that  this  choice  of  parameters  has  the  following  special  features,  which 
are  not  common  to  the  best  operating  thermometer  devices:  very  low  drive,  or 

carrier  amplitude;  drive  frequency  near  the  fourth  harmonic  of  the  Josephson 
relaxation  oscillation  frequency,  which  appears  to  lead  to  phase  locking  between 
the  two  oscillations,  and  sharp  voltage ■ "spiking"  associated  with  the 
negative-resistance  intervals  when  (l+ycos<j))  is  negative. 

Thus,  a value  for  oj  of  45.68  for  this  choice  of  R and  Ic  would  correspond  to 

45.68  x 4.836  = 220.9  kHz 

the  observed  average  relaxation  frequency  is  101/421  of  that  value;  this 
corresponds  roughly  to  the  ratio  between  a =11  and  cj  = 45.68. 

There  is,  however,  a perturbation  associated  with  the  amplitude  factor  k for 
the  drive  term.  The  theory  by  Peterson,  Soulen,  and  Giffard  gives  a value  for 
this  correction  to  the  predicted  relaxation  frequency 

(1  - J0* 2/(ll)2)1/2 

where  the  argument  of  the  Bessel  function  J0(k)  is  1 in  the  case  we  have 
chosen.  The  corrected,  predicted  relaxation  frequency  would  then  be 
53.0651  kHz;  we  find  in  simulation  the  value  52.9970.  The  reason  for  this 
discrepancy  of  0.13%  is  not  clear  — the  observed  value  is  lower  than  the 
predicted. 

Results  of  the  Integration 

When  started  from  <}>  =0,  dtp  /d  t =0  , a point  which  is  off  the  stable 

trajectory,  the  solution  evolves  quasi-periodically  for  about  1000  drive 
cycles,  then  settles  into  the  stable  periodicity,  reproducing  itself  to  about 

2 parts  in  10  , in  both  $ and  dtp  / dt • It  is  not  clear  whether  failure  to 

cycle  more  closely  than  this  represents  a fluctuation  associated  with  phase 
locking,  in  the  physical  system.  Both  exact  and  approximated  Jacobian 
routines  have  been  tried  at  varying  levels  of  precision,  with  similar  results. 
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We  find  the  numerical  solution,  in  double  precision,  to  be  multiply  periodic, 
with  a prominent  subharmonic  at  421  times  the  drive  period  2tt/w  , and  at  101 
times  the  relaxation  period.  This  determines  the  average  relaxation 
frequency,  from  the  repetition  rate,  to  be  101/421  x 45.6800  very  precisely 
(to  the  order  of  parts  in  10  ).  There  is  a family  of  nearly-periodic 
structures  appearing  ; "quasi-subharmonics"( ? ) reflected  in  the  patterns  of 
modulation  of  the  high-frequency  spikes  in  d$/dt. 

It  is  not  clear  at  this  early  stage  of  the  analysis,  whether  this  rational 
subharmonic  relation  is  a consequence  of  frequency-locking  of  the  relaxation 
oscillation,  via  subharmonics,  to  the  nearest  subharmonic  of  the  drive 
frequency.  We  note  that  for  the  parameters  chosen,  the  drive  frequency  is  a 
little  more  than  four  times  the  relaxation  frequency,  and  hence  " band 
separation"  is  incomplete.  Also,  the  observed  relaxation  frequency  is  offset 
towards  a lower  value  than  that  suggested  by  the  simple  theory,  e.g.,  there 
may  be  a "pulling"  of  the  relaxation  frequency  by  the  frequency  locking. 

Operation  of  the  program  at  single  precision,  in  other  parameter  ranges, 
indicates  that  the  spikes  originate  with  the  increase  of  y to  about  one;  when 
y reaches  a value  about  1,  at  this  drive  amplitude,  the  spikes  may  be  negative 
in  sign  as  well  as  positive;  at  1.5  they  can  occur  in  multiples  as  well  as 
singly. 

Fig.  23  $ , the  superconducting  phase-difference  across  the  tunneling  device, 
is  shown  modulo  2tt  as  a function  of  time.  Note  that  on  the  up-going 
transitions  between  steps,  in  the  relaxation  oscillation,  only  a relatively  few 
points  are  caught  in  the  § transition  region  near  tt  . The  function  is  a 
continuous  staircase  moving  upward  off  the  plot. 

Fig.  24  $ , the  instantaneous  voltage,  is  shown  on  a very  much  reduced 

scale.  There  are  sharp  positive-and  negative-going  excursions  (to  positive  800 
and  to  negative  600)  during  the  very  brief  time  intervals  when  $ is  near  tt  . 
Other  plots  indicate  that  its  mean  value  lies  around  4 to  5 on  the  stable 
manifold  points.  Evidently  the  spikes  contribute  a substantial  part  of  the 
mean  voltage  ~ 10  or  11. 

Fig.  25  sin<£  , the  instantaneous  current  in  the  junction  loop,  is  displayed 

for  the  same  points.  Roughly,  it  appears  as  a sinusoid  at  the  drive 
frequency,  added  to  a linear  ramp  from  the  relaxation  oscillation.  (Note  that 
this  ramp  is  not  derived  from  the  numerical  "chopping",  modulo  2tt  , of  Fig.l, 
and  is  smoothly  continuous  across  those  drops);  the  drops  in  Fig.  25  represent 
quite  drastic  reversals  of  the  supercurrent  direction  through  the  junction, 
including  some  "extra"  reversals  going  to  the  current  limits,  associated  with 
the  large  spikes  in  the  plot  of  the  instantaneous  voltage.  In  the  case  of  the 
Josephson  point  contact  facing  into  an  inductance,  the  phenomena  represent 
pulses  of  supercurrent  sharpened  by  the  inductance.  Closer  examination 
reveals  that  Fig.  25  actually  "wraps"  Fig.  23  around  a cylinder,  with  the 

"jump"  regions  of  very  rapid  current  reversal  occurring  on  the  " far"  side,  as 

viewed  in  projection. 
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Figs.  26  and  27  z = 6 <j>  + <j>  + y sin  <j>  the  coordinate  which  was 

suggested  as  an  alternative  separation  variable,  is  shown  plotted  against  4>  , 
for  2 sets  of  points  taken  at  time  intervals  corresponding  to  tt/w  and  0.1  tt  / u> 
from  the  time  zero  of  the  sin  cot  driving  term,  respectively.  (a  modified  type 
of  Poincare  plot).  Two  branches  of  the  slow  manifold  are  seen;  the  points 
corresponding  to  tt/co  wobble  upward  on  the  left-hand  branch  to  the  maximum,  then 
jump  across  the  central  "unstable"  zone.  The  points  corresponding  to  0.1  tt/(o 

....  wobble  downward  on  the  right-hand  branch,  then  jump  across  to  the  left 

[with  a lower  minimum  slope  than  the  up-going  jump  points]. 

Figs.  28,  29,  30  and  31  display  the  same  data  as  a function  of  sin  <j>  . We  note 

that  <p  appears  to  be  a cyclic  function  of  z,  modulo  2u ; however  we  can  only 

conclude  from  this  that  it  is  a continuous  function  of  z. 

Figs.  28  and  29  are  representative  of  fully  stabilized  cycles  in  the  solution. 
Figs.  30,  31,  and  32  are  data  which  include  the  early  transient  behavior  as  the 
system  is  started  from  (0,0)  at  t=0. 

Fig.  33  which  shows  z as  a function  of  t,  should  be  compared  against  Fig.  23 
and  25.  It  will  be  noted  that  the  drops  in  Fig.  33  are  only  those  associated 
with  the  arithmetic  operation  of  reducing  <j>  to  values  modulo  2 tt  . z in  this 
plot  is  simply  the  sum  of  terms  phi,  phidot,  and  sin  <j>  , evaluated  from  the 

returned  <j>  and  cj>  of  the  integrator  program,  with  <j>  expressed  modulo  2 tt  . 

z as  a function  of  t evolves  smoothly  as  a linear  ramp,  with  an  added  sinusoid 
at  the  "carrier"  frequency.  The  slope  of  the  ramp  reflects  the  relaxation 

frequency.  This  smooth  evolution  takes  place  through  the  " jump" 

discontinuities  in  the  time  derivatives  of  the  other  functions  plotted;  its 
behavior  confirms  the  general  internal  consistency  of  those  solutions. 

Thus  z vs  t is  like  Fig.  23  in  portraying  the  continous  evolution  of  a 
"rotating"  solution;  however,  its  rate  is  uniform,  and  does  not  exhibit  the 
cycle  of  the  relaxation  oscillation. 

Long-Time  Behavior 

Figs.  34  and  35  illustrate  Long-time  or  Very  Slow  Behavior 

If  the  results  illustrated  in  Fig.  3 and  Fig.  5 are  extended  by  plotting 
only  one  point  per  drive  cycle,  and  greatly  compressing  the  time  scale  (3000 
drive  cycles,  or  about  12  milliseconds,  is  the  duration  of  the  plot  shown  here) 
the  result  is  an  interesting  "stroboscopic  imaging"  of  the  long-term  motion. 
After  the  first  1000  drive  cycles,  the  motion  is  fully  periodic  and  the  diagram 
repeats  itself  to  parts  in  10  . 

The  upward  "return  tracks"  record  the  intervals  at  which  extreme  values  of 
voltage  spikes  recur  (the  subharmonic  421/101). 
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Fig.  35  is  the  same  type  of  plot,  but  calculated  with  y = 0.5  instead  of  1.5. 
The  "stroboscopic  path"  is  greatly  smoothed,  and  more  nearly  sinusoidal.  In 
this  case,  voltage  "spiking"  is  not  present.  The  maximum  value  of  phidot 
encountered  in  this  run  is  about  14. 

Careful  examination  of  this  diagram  reveals  that  the  pattern  at  the  right  hand 
end  is  starting  to  repeat,  after  time  = 330.4  on  the  time  scale  of  the  plot. 
Examination  of  printout  shows  that  this  doubly-periodic  solution  is 
reproducing  itself  to  within  3 or  4 parts  in  10  . 

The  repetition  period  corresponds  to  2402  drive  cycles,  or  577  periods  of  the 
relaxation  oscillation.  The  high  rational  "winding  ratio"  may  be  used  to 
estimate  the  numerical  value  of  the  relaxation  period,  averaged  over  the  long 
repetition  period  and  expressed  in  units  of  2t:/w 


Trel  = 2402  = 4.16291 

577 

For  comparison,  the  "autonomous"  relaxation  period  expressed  in  the  same  units 

w = 45.680  = 4.1527273 

11.0 

when  corrected  by  the  Gif f ord-Soulen-Peterson-Sanders  factor 

a 

= ( .9975775)'1 

/ a^-J  0Z (k) 

for  a = 11  k = 1 yields  the  value  4.16281. 

We  conclude  from  this  that  the  "autonomous"  relaxation  period,  corrected  by 
this  Bessel-function  formula,  agrees  accurately  with  the  average  relaxation 
period  seen  in  the  simulated  solution  of  the  second-order  equation,  for  y =0.5 


Figs.  36  and  37  show  confirmation  between  a result  from  Sanders,  on  a problem  of 
critical  variation,  and  the  same  answer  as  calculated  by  this  program.  For  this 
test  the  frequency  w is  set  at  45.68128  (slightly  different  from  the  values  used 
in  the  preceding  series  of  figures),  and  the  initial  value  of  $ is  varied  in 
steps  of  less  than  2 parts  in  10^. 

The  initial  value  of  phi  is  chosen  as  follows: 

a)  6.374423 

b)  6.374424 

c)  6.374425 

d)  6.374426 

e)  6.374427 

f)  6.374428 

This  variation  produces  changes  of  the  order  of  1 in  the  output  waveform  (phi 
versus  time),  which  is  that  corresponding  to  a triple  crossing  of  the  region 
where  1 + y cos  <p  < o. 
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In  Fig.  37  the  horizontal  lines  represent  the  values  of  <j>  for  which  1 + y cos<{) 
vanishes . 

Summary  of  Observations 


These  tests  of  quiet  and  "jump"  cases  have  produced  a variety  of  graphical  and 
numerical  output.  We  summarize  here  a few  of  the  more  interesting  results. 

Quiet  cases 

o For  y < 1 ( at  least  as  high  as  0.8)  all  waveforms  are  smooth, 

regardless  of  the  amplitude  of  the  drive  signal.  There  is  a gradual 
transition  resembling  that  from  amplitude-to  frequency-modulation,  as 
the  drive  amplitude  is  increased. 

o A useful  memory  device  is  to  think  of  <f>  as  an  argument  function 
confined  to  move  on  a cylinder.  Sin  <j>  , the  instantaneous 
supercurrent,  is  then  visualized  as  a projection  on  a plane  parallel  to 
the  axis  of  the  cylinder. 

o A type  of  mathematical  "demodulation”  or  separation  of  the 

low-frequency  relaxation  oscillation  from  the  high-frequency  "carrier" 
is  easily  achieved  by  making  the  return  interval  of  the  integrator 
program  coincide  with  the  period  of  the  drive  signal. 

o Tests  at  a low  value  of  (a-1 ) give  good  agreement  with  the 

Bessel-function  dependence  of  relaxation  frequency  (period)  on  drive 
amplitude  coefficient  k . 

"Jump"  case 


o As  expected  from  first-order  theory,  the  simulation  exhibits  "jumps" 
between  parts  of  the  slow  manifold.  These  rapid  transitions  can  be 
accurately  traced  out  in  full  detail,  provided  by  the  second-order 
equation.  They  are  accompanied  by  sharp  voltage  pulses. 

o The  solution  is  observed  to  be  doubly  periodic  at  high  precision,  with 
the  relaxation  oscillation  at  the  421/101  subharmonic  of  the  drive 
frequency.  Since  the  drive  signal  has  been  selected  near  the  4th 
harmonic  of  the  relaxation  frequency,  we  suspect  that  phase-locking  is 
occurring  because  of  the  strong  overlap  between  the  "audio"  and 
"carrier"  bands.  Thus  it  may  be  difficult  to  compare  the  observed 
repetition  period  of  the  relaxation,  with  a theoretical  prediction 
based  solely  upon  values  of  a , k , and  y . Unlike  the  quiet  cases  which 
were  considered  above,  the  solution  may  have  a relaxation  frequency 
which  wanders  from  a constant  value  during  the  "long"  421  repetition 
period,  while  averaging  to  a mean  value  which  is  frequency-locked  to 
the  drive  frequency. 

o Plots  of  z versus  phi  exhibit  sharp  details  of  the  "up-going"  and 

"down-going"  transitions  between  parts  of  the  slow  manifold.  If  phi  is 
shown  in  cylindrical  projection,  the  connection  with  current  reversals 
going  to  the  critical  limits  can  be  made  clear. 
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o Plots  of  z versus  time  exhibit  the  uniform  advance  of  a smooth  function, 
apparently  sinusoidal  at  the  drive  frequency,  atop  a highly  linear  ramp; 
this  in  spite  of  the  fact  that  its  three  constituent  terms  show 
discontinuities  in  their  derivatives.  This  suggests  the  possibility  of 
further  simplifications,  in  pursuit  of  a numerically  simple  routine  for 
evaluating  relaxation  period  (frequency). 

Conclusions  as  to  the  Accuracy  of  the  Integrator 

o Figs.  36  and  37  demonstrate  the  independence  of  the  solution  method  from 
details  of  algorithm,  programming,  and  machine;  and  also  the  high 
sensitivity  of  both  NBS's  and  Sanders’  algorithms  to  critical  changes  in 
initial  conditions.  The  figures  are  an  illustration  of  how  small  changes 
in  the  initial  conditions  can  generate  dramatic  changes  in  the  behavior  of 
the  wave  form  as  <j)(t)  traverses  the  region  where  1+y  cos  $ < 0.  This  is 
an  instance  of  the  changes  from  ’dips'  to  'slices'  which  were  studied  by 
Littlewood  for  the  case  of  the  van  der  Pol  equation. 

A tolerance  of  10” is  needed  to  produce  Figs.  36  and  37.  Accumulation 
of  error  is  a well  known  phenomenon  in  numerical  solution  of  o.d.e.'s. 

The  accumulation  of  error  eventually  gives  a completely  wrong  answer.  An 
accumulated  error  of  one  part  in  10^  would  alter  the  traversal  of  the  1+y 
cos  cp  < 0 region  (as  can  be  seen  in  Figures  36  and  37). 

o The  system  is  not  "stiff".  Considerations  of  accuracy  restrict  the  step 
size.  We  obtained  the  most  efficient  integrations  by  using  the  Adams 
Moulton  integrator  coupled  with  Newton's  method  for  the  nonlinear 
equations.  However,  given  sufficient  machine  precision  any  competent 
modern  integration  will  correctly  follow  these  solutions.  This  was 
confirmed  by  Sanders,  et  al , who  obtained  some  similar  results  using 
entirely  different  software. 

o The  observed  doubly-periodic  cycling  of  the  "jump"  case  for  gamma  =1.5  is 
followed  to  parts  in  10^  for  the  observed  output  values  of  phi  and  phidot, 
after  an  initial  "settling"  period  of  about  1000  drive  cycles.  We 
associate  this  "settling"  period  with  the  approach  of  an  "off-manifold" 
solution  to  the  asymptotically  stable  form. 

o The  results  just  mentioned  were  obtained  only  at  double-precision  10”^ 
tolerance  requirement.  Improving  the  tolerance  to  parts  in  10-*-0  did  not 
produce  distinguishable  improvements.  By  contrast,  single  precision 
(parts  in  10^)  produced  "time-slippage"  which  threw  off  the  accuracy  of 
the  long-time-period  cycling.  Although  general  features  of  the 
single-precision  runs  were  similar  to  the  double-precision  ones,  they 
clearly  lacked  the  full  detail  and  accuracy.  (We  note  that  errors  in  phi 
and  phidot  are  cumulative  with  time,  so  that  an  error  of  one  part  in  10^ 
per  cycle  can  accumulate  to  one  part  in  10^  after  1000  cycles). 

These  results  are  in  accord  with  expectations,  since  the  double-precision 
integrator  is  expected  to  perform  distinctly  better  in  the  rapid 
transitions  between  the  "slow  manifold"  and  the  "jump"  intervals. 

For  the  "quiet"  cases  where  gamma  is  less  than  1,  the  distinction  in 
behavior  between  single-  and  double-precision  routines  is  expected  to  be 
less  sharp,  since  they  do  not  have  to  adjust  to  the  "jump”  transitions. 
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CONTINUING  DIRECTIONS  OF  THE  INVESTIGATION 


o With  the  reliability  of  the  integrator  now  confirmed  by  severe  tests,  we 
wish  to  apply  it  to  some  major  unsolved  problems  of  the  experiments.  One 
such  is  the  behavior  of 


<<£>  = lim  _1  /T  <£(t)dt 

X-xx>  'p  O 

vs  k for  a driven  y>l  case.  Sanders  t11]  gives  an  estimate  of  the  period 
T0+AT  for  a y>l  autonomous  case,  and  the  formula 

2tt  w 

T = >V-  JQZ  (k) 
for  a y<l  driven  case. 

We  are  trying  to  learn  more  about  the  observed  anomalous  "doubling"  in 
the  driven  y>l  case.  We  had  to  have  an  utterly  reliable  integrator 
before  attempting  this. 

Foremost  in  interest  will  be  to  run  with  a wide  range  of  drive  amplitude 
values,  up  to  k = 8 or  9,  for  values  of  gamma  between  1 and  2;  the  output 
quantity  of  interest  will  be  phi  divided  by  t-t0,  evaluated  to  a part  in 
10^,  for  comparability  with  the  experimental  observations,  which  have 
shown  examples  of  "doubling”  at  the  parts-in-lO^  level,  not  in  accord 
with  a simple  Bessel  expression. 

o Return-mapping  techniques  will  be  tried,  in  order  to  examine  questions  of 
frequency  stability,  especially  the  possibility  of  chaotic  behavior. 
Chaotic,  or  erratic  behavior  of  the  observed  relaxation  frequency  would 
be  expected  to  interfere  with  the  desired  high-fidelity  operation  of  the 
Josephson  thermometer  as  a voltage-to  frequency  converter  of  thermal 
(Nyquist)  noise. 

o An  asymptotic  expansion  formula  derived  by  Sanders  L11]  may  provide  a key 
to  developing  Fourier  analysis  of  the  "quiet”  waveforms,  such  as  Figs. 
16-21.  This  will  be  explored  in  further  detail. 

o The  programs  will  be  rendered  transportable  (with  consideration  being 
given  to  its  possible  future  use  on  smaller  machines  than  the  Univac 
1100/82).  Jill  Novotny  is  currently  installing  it  on  the  Catholic 
University  DEC  10. 
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APPENDIX  B 


CIRCUIT  ALGEBRA 


TANK 

CIRCUIT 


from  dc 
and  rf  feed 


1 


'rO 

■GSW — 


M 


RSQUID 

THERMOMETER 

LOOP 


dc 


Point  Contact 

(parasitic)  capacitance  of  this 

Resistance  associated  with  the  quasiparticle  conduction 
Single  coupling  loop  (inductance  at  nanohenry  level) 

Shunt  Resistor,  typically  of  silicon  bronze 
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Josephson  relations  in  the  point-contact  branch  (an  active  device) 


Vj  = 1/2it  • h/2e  • d^/dt  instantaneous  Josephson  emf  (la) 

Isc  = Ic  s;*-n  ‘Kt)  = tunneling  supercurrent  (lb) 

Current  division: 

Ij[n  3 Ij_  + Ie  at  input  terminal  to  RSQUID  loop  (2) 

II  = Ij  + Ic  + I2  at  input  to  point  contact  (3) 

Displacement  current  in  capacitance 

Ic  = C • dVj  /dt  (4) 

Voltage  closure  around  the  emf’s  of  the  main  RSQUID  loop 

Vj  = M dlrf/dt  + L dl1/dt  + R^Ie  (5) 


We  assume  that  the  tank  circuit  functions  as  a constant-current  rf  source, 
and  that  the  dc  bias  feed  comes  from  a constant-current  source,  so  that 

Iin  = nlrf  + Idc  n < 1 , where  Irf  = irf  sin  ut  (6) 

Here  Irf  is  the  rf  current  resonating  in  the  tank  circuit,  and  n 
represents  the  fraction  of  this  which  is  diverted  into  Ij_n 


We  substitute  into  (5)  and  combine  these  relations  into  a single  equation  for 
the  emf’s  around  the  main  RSQUID  loop.  With  regrouping  of  terms,  we  obtain 
the  full  equation  for  <j>  as  a function  of  time. 


1/ 2tt  • h/2e  { - L Cj  d<{>3/dt3  + [ RgCj  - L/Rj  J d$2 / dt2 

+ [ 1 + Re/Rj  - Ic  L cos  9]  d(j>/dt  } + Re-^-c  s:i-n  9 
R^  Idc  + [ M a)  + n Re  j irf  sin  ut  (7) 
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The  third  derivative  term,  and  the  second  term  within  the  second  bracket  are 
neglected.  We  take  four  new  definitions  of  parameters: 


0 E Re  Cj  - L/Rj 

Y = - Ic  L 


a = (Rg  Idc  )/(  Re  Ic  ) = Idc  / Ic  dc  bias 

[M  + n Re]  irf 

kcu  = rf  drive  amplitude 

Re  Ic 


and  note  that  there  is  no  loss  of  generality  in  assuming  y positive,  except 
that  we  must  recall  this  reversal  of  sign  when  re-interpreting  emf’s  on  the 
inductance  L.  This  achieves  the  expected  positions  of  the  coefficients. 

If  in  addition  we  scale  the  frequency  m to  the  plasma  frequency  of  the 
RSOUID  loop,  which  is  set  by  the  coefficient  of  the  sine  term  (the  product  of 
the  critical  current  Ic  and  the  shunt  resistor  R€): 

2tt  ojc  = 2e/h  • Re  Ic  =1 

we  obtain  the  greatly  simplified  general  form  for  the  RSJ  equation  of  the 
RSQUID  loop. 


3<J>  + (1+ycos  !p)(f>  + sin  <j)  = a + ka)  sin  cot  (8) 

Here  we  have  tacitly  assumed  that  Re/Rj  <<  1 , so  that  in  fact  we 

neglect  entirely  the  quasiparticle  tunneling  current  through  the  point  contact 
at  the  very  low  temperature  of  the  RSQUID  thermometer.  This  enables  us  to 
drop  the  subscript  on  Rg  and  to  write  it  as  R in  the  scaling  description 
starting  on  p.  11. 

When,  in  addition,  we  assume  a positive  value  of  .001  for  0 , we  are 
implicitly  forcing  a "smoothing"  time  constant  upon  the  system  response  during 
the  short  "jump"  intervals: 

TRC  = Re  C ~ 0 • 2tt/coc  = .21  microseconds 


This  detailed  circuit  model  was  derived  in  a slightly  different  manner  by  D. 
van  Vechten  (private  communication,  1980)  and  has  been  partially  published  in 
references  (7),  (8),  and  (9). 
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